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Abstract 

A genetic algorithm approach suitable for solving multi-objective optimization problems is described and 
evaluated using a series of aerodynamic shape optimization problems. Several new features including 
two variations of a binning selection algorithm and a gene-space transformation procedure are included. 
The genetic algorithm is suitable for finding Pareto optimal solutions in search spaces that are defined by 
any number of genes and that contain any number of local extrema. A new masking array capability is 
included allowing any gene or gene subset to be eliminated as decision variables from the design space. 
This allows determination of the effect of a single gene or gene subset on the Pareto optimal solution. 
Results indicate that the genetic algorithm optimization approach is flexible in application and reliable. 
The binning selection algorithms generally provide Pareto front quality enhancements and moderate 
convergence efficiency improvements for most of the problems solved. 
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Nomenclature 

set of scalar objective functions 
optimal set of F values 
scalar objective function 
maximum value of f 
n th GA generation (active file) 
integer ranking array 

seed value used in the random number generator R(0,1) 

user-specified integer parameter controlling which selection algorithm is used 

user-specified integer parameter controlling the gene-space transformation option 

fixed number of chromosomes in each GA generation 

number of genes in each chromosome 

number of scalar objective functions 

user-specified vector with four elements that controls modification operator usage 
user-specified parameter controlling probability that a specific gene will be modified using 
perturbation mutation operator (0 < p., < 1) 

user-specified parameter controlling probability that a specific gene will be modified using 
global mutation operator (0 < p 2 s 1) 
basis vector matrix with elements r,j 

random number generator that returns a random value between 0 and 1 
unitary transformation matrix with elements u, j 

/ th gene from the y' th chromosome and the n th GA generation 
y' th chromosome from the n th GA generation 
user-specified maximum limit on the /' th gene 
user-specified minimum limit on the / th gene 

user-specified parameter controlling the size of the perturbation mutations (0 < /3 < 1) 
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subscripts 

/ 

j 

k 


gene or decision variable index 
chromosome index 
objective function index 


superscripts 

n GA generation or population index 

t temporary chromosome values obtained after selection but before modification 


Background 

Numerical methods for optimizing the performance of engineering problems have been studied for many 
years. Perhaps the most widely used general approach involves the computation of sensitivity gradients. 
These methods— called gradient methods — have been utilized to produce optimal engineering 
performance in a wide variety of different forms. The reliability and success of gradient methods generally 
requires a smooth design space and the existence of only a single global extremum or an initial guess 
close enough to the global extremum that will ensure proper convergence. 

In contrast to gradient-based methods, design space search methods such as genetic algorithms (GA) — 
also referred to as evolutionary algorithms (EA)— offer an alternative approach with several attractive 
features. The basic idea associated with the GA approach is to search for optimal solutions using an 
analogy to the theory of evolution. The problem to be optimized is parameterized into a set of decision 
variables or genes. Each set of genes that fully defines one design is called an individual or a 
chromosome. A set of chromosomes is called a population or a generation. Each complete design or 
chromosome is evaluated using a “biological-like” fitness function that determines survivability of that 
particular chromosome. For example, in aerospace applications, the genes may be a series of geometric 
parameters associated with an aerospace vehicle that is to be optimized for payload delivered to orbit, 
aerodynamic performance, or structural weight. The fitness function takes as input all the geometric 
parameters and returns a numerical value for the fitness— the size of the payload, the aerodynamic 
performance, or the structural weight. 

During solution advance (or “evolution” using GA terminology) each chromosome is ranked according to 
its fitness vector— one fitness value for each objective. The higher-ranking chromosomes are selected to 
continue to the next generation while the lower-ranking chromosomes are not selected at all. The newly 
selected chromosomes in the next generation are manipulated using various operators (combination, 
crossover, or mutation) to create the final set of chromosomes for the new generation. These 
chromosomes are then evaluated for fitness and the process continues— iterating from generation to 
generation — until a suitable level of convergence is obtained or until a specified number of generations 
has been completed. 

Constraints can easily be included in the GA optimization approach either by direct inclusion into the 
objective function definition— a so-called penalty constraint— or, more commonly, by including one or 
more constraints into the fitness function evaluation procedure. For example, if a design violates a 
constraint, its fitness is set to zero, that is, it does not survive to the next generation. Because GA 
optimization is not a gradient-based optimization technique, it does not need sensitivity derivatives. It 
theoretically works well in non-smooth design spaces containing several or perhaps many local extrema. 

General GA details including descriptions of basic genetic algorithm concepts can be found in Goldberg , 1 
Davis , 2 and Beasley, et al . 3,4 Additional useful studies, which survey recent activities in the area of GA 
research including the presentation of model problems useful for evaluating GA performance, are given 
in Deb , 5 Van Veldhuizen and Lamont, 6 and Jimenez, et al . 7 
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A disadvantage of the GA approach is expense. In general, the number of function evaluations required 
for a GA optimization process to converge, exceeds the number of function evaluations required by a 
finite-difference-based gradient optimization (see the results presented in Obayashi and Tsukahara 8 and 
Bock 9 ). This situation is offset, to an extent, by the ease with which GAs can be implemented in parallel 
or distributed computing environments. 

Despite being relatively new, genetic algorithms have already been applied in many applications over a 
broad range of engineering fields. A brief survey of single discipline applications is presented in Holst 
and Pulliam , 10 and will not be discussed further in this report. 

Other applications involving GA search methods have been made in the area of multi-objective or multi- 
discipline optimization, that is, optimization problems in which two or more objectives are simultaneously 
and independently optimized. These methods, referred to as MOGA (multi-objective genetic algorithm) 
methods, are especially attractive because they offer the ability to directly compute so-called “Pareto 
optimal sets” in a single computation instead of the limited single design point traditionally provided by 
other methods. 

The Pareto optimal set or Pareto front, as it is commonly called, includes optimal solutions for each of 
the individual objectives, as well as a range of tradeoff solutions in between, which are themselves 
optimal solutions. Providing a range of solutions to a multi-objective optimization problem is a powerful 
approach because it allows the designer to see the effect of decision variable variation on the design 
space in the form of optimal tradeoffs. Thus, the designer can choose individual objective weighting 
factors after their full influence is quantitatively known. 

In recognition of the importance of the MOGA approach, many theoretical developments have been 
published in recent years. In particular, Van Veldhuizen and Lamont 6 and Zitzler, et al . 11 present 
reasonably complete surveys of MOGA methodology with a special emphasis on how to compare GA 
performance from one algorithm variation to another. Other researchers, including Deb , 5 Deb, et al ., 12 
Van Veldhuizen , 13 Lohn, et al . 14 and Holst and Pulliam , 15 have developed and/or utilize a suite of test 
problems as a standard in evaluating MOGA convergence efficiency as well as the accuracy of the final 
Pareto front. A key aspect of Pareto front development is diversity. Does a particular MOGA produce 
good coverage over the entire Pareto front or are some regions poorly resolved while other regions have 
high levels of undesirable clustering? This issue is studied in Laumanns, et al . 16 and Deb, et al . 17 

Still other MOGA issues are associated with archive strategies. As the Pareto front develops, many 
solutions are found that lie on the Pareto front that cannot be retained within the fixed population size of 
many schemes. How to retain or archive this information for the benefit of the MOGA while maintaining 
reasonable bounds on the archive file has been studied in Knowles and Crone . 18,19 One last example of 
MOGA research lies in the area of an interactive algorithm development. Parmee, et al . 20 present a 
MOGA approach that allows changes to be made in GA operators as well as in objective definitions as 
the MOGA advances from generation to generation. 

One area of MOGA research that is even more voluminous than the area of theoretical developments is 
that associated with GA multi-objective applications. Examples of multi-objective optimization 
applications include airfoil optimization by Marco, et al ., 21 Naujoks, et al ., 22 Quagliarella and Della 
Cioppa , 23 Vicini and Quagliarella , 24 Hamalainen, et al . 25 and Epstein and Peigin , 26 missile aerodynamic 
shape optimization by Anderson, et al ., 27 and wing optimization by Anderson and Gebert , 28 Sasaki, et 
al ., 29 Oyama , 30 Obayashi, et al ., 31 and Ng, et al . 32 

Additional MOGA examples in the area of turbomachinery optimization include rocket engine turbopump 
design by Oyama and Liou 33 and compressor blade design by Benni 34 and Oyama and Liou . 35 In some of 
these examples the multiple objectives were obtained by considering two different aerodynamic design 
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points. In others the multiple objectives involved different disciplines including aerodynamics, structures, 
controls and/or electromagnetics. 

One last area of GA research that bears mention is in the area of hybrid methods, the utilization of GA 
optimization in conjunction with another type of optimizer. This is an attractive area of research because 
GA methods are particularly good at finding a global extrema, but not particularly good in converging the 
optimal solution to tight levels of accuracy. Briefly, several examples where hybrid approaches have 
been utilized include Rai 36 and Madavan 37 for single objective problems and Giotis, et al., 38 Tursi, 39 Vicini 
and Quagliarella 40 and Brown and Smith 41 for multi-objective problems. 

Various definitions and the multi-objective genetic algorithm used in the present study are described next. 
Details associated with each of the operators, including selection, passthrough, random average 
crossover, perturbation mutation, and mutation are presented. MOGA effectiveness is then evaluated 
using several three-dimensional aerodynamic shape optimization problems. 


Problem Statement: Single Objective Optimization 

A single-objective optimization problem can be stated as follows: Let f be a scalar function of N G 
independent variables, x, , defined on some domain Q 


f = f(X) = f{x^...,x h ...,x NG ) (la) 

In this notation X is the vector of design space decision variables. The maximum value of f , indicated by 
f* , is obtained by finding the values of X =X* such that* 

r = max{f} = f(X*) = f{xl...,x*,...,x* Ne ) (1b) 

The above maximization operation is subject to N co conditions or constraints indicated by 

c n (X)<0 n = 1,2,...,A/ co (1c) 

The constraints placed on the decision variable vector X by Eqs. (1c) essentially serve to limit the design 
space within Q for which the optimal solution is sought. 


Problem Statement: Multi-Objective Optimization 

For optimization problems involving more than one objective, which are simultaneously and independently 
optimized, the situation is more difficult. This is because each objective must play a role in determining 
the optimal solution. In the optimization process, conflicts might arise among the various objective 
functions, that is, the optimal values of each individual objective, in general, will not occur for the same 
decision variable vector, X . As a result, the “optimal solution” for a multi-objective optimization problem is 
typically a range or a set of solutions that represent tradeoffs in objective space. 

To determine when one solution is better than another for multi-objective optimization problems the 
concept of dominance is utilized. 1 A vector U = U (u- i ,...,u i ,...,u N ) is said to dominate another vector 
V = V(v 1 ,...,v, ,...,v w ) if and only if u , >v ( for all / and there exists at least one value of / such that 


* For the purpose of simplifying the discussion of algorithmic details, maximization is generally assumed. 
The logic for minimization is a straightforward modification and will not be discussed. 
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Uj >Vj. A vector defined on some domain £2 that is not dominated by any other vector defined on £2 is 
said to be non-dominated on £2 . 

A multi-objective optimization problem can be stated as follows: Let F be a set of N a scalar objective 
functions, f k , each dependent upon the same decision variable vector, X, which is defined on some 
design space £2 


F = F(f,(X),...,f fe (X) f Ng (X)) (2a) 

As above, the decision variable vector X consists of N G independent components. The multi-objective 
optimization problem involves finding the set of X = X* that produces non-dominated values for F = F* on 
£2 . This set of values F* is called the Pareto optimal set or the Pareto front. 

For each F the constraints 


c n (X)<0 n = 1,2 N co (2b) 

must be satisfied. Existence of these constraints serves to limit or reduce the size of £2 for which the 
optimal solution is sought. 


Genetic Algorithm 


The genetic algorithm optimization procedure utilized to solve the multi-objective optimization problem, as 
described by Eqs. (2), is now presented. It is closely related to the MOGA optimization procedure 
presented in Holst and Pulliam. 15 As mentioned in the introduction, the general idea behind GA 
optimization is to discretely describe the design space using a number of decision variables, x,. In GA 
parlance these parameters are called genes, and the / subscript refers to the gene number. Each set of 
genes that leads to the complete specification of an individual design, that is, each decision variable 
vector, X , is called a chromosome and is indicated by 


*7 


:X"(x" y ,...,xT,...,x£ sJ ) 


( 3 ) 


where the j subscript, added to X, identifies the chromosome number. In addition, the j subscript has 
been added to each gene, to indicate its chromosome membership. The n superscript has been added to 
indicate the GA generation number, which is iteratively advanced as the solution converges. Thus,X y is 
the y th chromosome for the n th generation that consists of N G genes. 


For aerodynamic shape optimization problems, the design space genes are typically a series of geometric 
parameters, for example, airfoil thickness and camber and/or wing sweep, twist and taper. For many GA 
applications genes are computationally represented using bit strings and the operators used to 
manipulate them are designed to accommodate bit string data. In the present approach, following the 
arguments of Oyama, 42 Houck, et al. 43 and Michalewicz, 44 real-number encoding is used to represent all 
genes. 

The key reason for using real number encoding is that it has been shown to be more efficient in terms of 
CPU time relative to binary encoded GAs. 44 In addition, real numbers are used for all genes in the present 
implementation because the present engineering application involves decision variables that are best 
described using real numbers, for example, the geometric parameters in aerodynamic shape optimization. 
Thus, using real number encoding eliminates the need for binary-to-real number conversions. 
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Initialization 


Once the design space has been defined in terms of a set of real-number genes, the next step is to form 
an initial generation, G°, represented by 


G°=(X°,...,X°,...,X° c ) (4) 

where N c is the total number of chromosomes. Each gene within each chromosome is assigned an initial 
numerical value using a process that randomly chooses numbers between fixed user-specified limits. For 
example, the / th gene in an arbitrary chromosome is initialized using 

Xj = /-?(0,1)(x max — -^min, ) + -*min, (5) 

where x max and x min are the upper and lower limits for the / th gene, respectively, and ft(0,1) is a 
random number generator that delivers an arbitrary numerical value between 0.0 and 1 .0. 

The random number generator used in the present study requires an integer input— a seed value. If the 
integer is positive, the next number in the current random number sequence is returned. If the integer is 
negative, the random number sequence is reset to a new sequence. Utilization of the same negative seed 
value will always reset the random number generator to the same random sequence. Each new solution 
begins by resetting the random number generator using a single call to R(0,1) with a negative seed 
value — ISEED. All other calls to R(0,1) during that solution use a positive seed value. Thus, a solution 
can be repeated by simply using the same initial seed value or rerun to determine statistical variation by 
using a different initial seed value. 

Fitness evaluation 

After a generation is established— either the initial generation or any succeeding generation in the 
evolutionary process— fitness values, F", are computed for each chromosome using a suitable function 
evaluation. This is analogous to the objective function evaluation in gradient methods and is represented 
using 


F, n =F(X?) ( 6 ) 

For example, for a multi-discipline optimization (MDO) problem involving the simultaneous maximization 
of two separate and distinct objectives, f, and f 2 , the fitness evaluation represented by Eq. (6) consists of 
the following 


f"-W) 

f 2 n =f2(Xj) 


where, for example, the first function might be the aerodynamic drag of an aerospace vehicle 
(constrained to fly at fixed lift) and the second function f 2 might be the structural mass of the same 
vehicle. Once all the genes in a specific chromosome have been specified, f, can be evaluated using a 
suitable CFD solver to obtain the drag, and f 2 can be evaluated using a suitable structural analysis 
routine to obtain the structural mass. In this case, of course, the optimization problem would be one of 
minimization. 
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Ranking 


The purpose of the ranking operation is to determine a set of integer values for the ranking array, IR . One 
integer value is required for each chromosome. The ranking algorithm is quite different for multi-objective 
optimization problems relative to single-objective problems. As such, both situations will be discussed in 
this section. The ranking array values— once determined— are then used in the GA selection process. 

Single-Objective Optimization— The GA ranking algorithm for single-objective optimization problems is 
quite simple. Whichever chromosome has the highest fitness is ranked number one (IR = 1), whichever 
has the second highest fitness is ranked number two (IR = 2), and so on. This ranking algorithm, more 
formally stated, is given by 


ic = 1 

if(fP<fO)ic~ic + 1 jj = \...,N c 
IR 1 = ic 


j =\-,N C 


(7) 


where j and jj are special counters that range over all N c chromosomes in the current population or 
generation level. 

Multi-Objective Optimization — For multi-objective optimization problems, the ranking procedure is more 
complex and utilizes the concept of dominance, which was defined previously. 

A chromosome with a fitness vector F that is not dominated by any other fitness vector within the design 
space is said to be a non-dominated chromosome. The optimal solution set F* or Pareto front includes all 
solutions with fitness vectors that are non-dominated and that satisfy the constraints given by Eqs. (2b). 
The numerical approximation to the Pareto front must involve a suitably complete set of discrete solutions 
so as to describe the optimal values of each individual objective— these are the Pareto front endpoints — 
as well as the non-dominated tradeoff solutions in between the optimal values associated with each of the 
individual objectives. 

With the concept of dominance in hand, the actual ranking process can now be presented. There are a 
number of ranking procedures available for use in multi-objective optimization. The one used in the 
present study is called Goldberg ranking. 1 After all of the fitness values have been determined, each 
chromosome is checked for dominance. Those chromosomes that are non-dominated are given a number 
one ranking {IR = 1) and then temporarily removed from the current generation. The remaining 
chromosomes that are non-dominated are given a number two ranking (IR = 2) and then temporarily 
removed. This process continues until each chromosome has an integer value for the ranking array, IR. 
In general, with this approach, the number of different integer values contained within the ranking array 
will be small, at least small in comparison to N c , as many chromosomes will be ranked near the top with 
a value of 1 , 2, or 3. 

The above procedure, by itself, represents a legitimate algorithm for the ranking operation. However, 
there is a refinement that provides a considerable enhancement to MOGA convergence. Before 
describing this refinement, two additional concepts must be defined— the active file and the accumulation 
file. 

The active file is simply a specific name used for the current generation of chromosomes, that is 

G n =(x^,...,xy,...,x^) ( 8 ) 
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is the n th generation active file for the GA iteration process. As the GA solution evolves, the active file is 
always fixed in size at N c chromosomes. The first N a elements in the active file— for the present 
implementation— always contain the chromosomes that possess the maximum fitness values for each of 
the N 0 objectives. 

The accumulation file is the list of all non-dominated chromosomes that have been found from all 
generations combined during the current MOGA iteration. As the MOGA solution evolves, the 
accumulation file typically grows in size with more chromosomes being added after each new generation. 
As new non-dominated chromosomes are added to the accumulation file, old chromosomes that lose their 
dominance must be deleted, thus ensuring that the accumulation file contains nothing but number-one 
ranked chromosomes. Because each chromosome stored in the n th generation accumulation file is non- 
dominated, it serves to define the Pareto front or, at least, the n th generation approximation to the Pareto 
front. The use of an accumulation file makes sense only when N 0 > 2 . 

The multi-objective ranking procedure described above utilizes only the chromosomes in the active file, 
that is, each chromosome in the active file is ranked relative to the other chromosomes in the active file. 
But this philosophy potentially wastes a plethora of information because not every number-one ranked 
chromosome for multi-objective optimization problems can be retained from one generation to the next in 
the active file. 

This is where the refinement in the ranking routine enters. Once each chromosome in the active file is 
ranked using the standard routine, an additional test is performed to see if any number-one ranked 
chromosomes in the active file are dominated by any of the chromosomes in the accumulation file. If this 
is the case, then the ranking number associated with the newly ranked chromosome in the active file is 
decremented by one. This ensures that the ranking routine produces number one rankings that are 
number one in a global sense. 

Selection 

After the ranking array is established, with or without the accumulation file option, the GA algorithm 
passes to the selection process to determine which chromosomes will continue on to the (n + 1) st 
generation and which will not. In the present algorithm implementation, four different selection operators 
have been coded: (1) a simple technique called greedy selection, (2) a traditional tournament selection 
algorithm, (3) a bin selection algorithm based on the computation of arc-length along the Pareto front, and 
(4) a second bin selection algorithm based on subdividing the design space into a matrix of boxes. 

The first two selection algorithms can be used for one or more objectives; the fourth can be used for two 
or more objectives, and the third inherently requires two objectives. The selection algorithm actually used 
is controlled by a user-input parameter called ISLECT, which will be defined shortly. 

Regardless of the selection algorithm chosen, the selection process picks established chromosomes — 
either from the n th generation active file, G n , or from the accumulation file— and then stores the results in 
a temporary holding array G f . The chromosomes stored in G f are then used by the subsequent 
modification operators to create G n+1 , which will be discussed shortly. For the present study, results are 
computed using only the first, third and fourth selection algorithms. Nevertheless, all four selection 
algorithms will be discussed for completeness. 
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Greedy selection— This selection algorithm is quite simple. It selects chromosomes from the n th 
generation active file, that is, from X" . It is implemented using the following: 

m = 1 

if (IRf < it) then 

Y 1 = Y n 

A /r? A y 
m = m + 1 

if(m > A/ c )stop 
endif 

where each selected chromosome X^, is placed in a temporary holding array indicated by 

G f = {X\,...,X t m ,...,X t Nc ) 

Note how the highest ranked individuals in the n th generation are selected multiple times, individuals with 
average ranking are selected a small number of times, and individuals with the lowest rankings are not 
selected at all. This biasing toward individuals with the highest rankings is a key element in any GA. This 
baseline variation of greedy selection is implemented by setting ISELECT = 1. 

Tournament selection— The tournament selection algorithm is quite simple and is used widely— one 
variation or another— in many GA optimization applications. The present variation is implemented as 
follows: Using a random process, three chromosomes are chosen from the n th generation active file, that 
is, from X". The ranking array values, !R, are then compared. The chromosome with the highest ranked 
array value is placed into a temporary holding array, G f . In case of ties the chromosome selected first is 
chosen. Next, three new chromosomes are chosen from X", and their ranking array values are 
compared. As before, the chromosome with the highest ranking is added to G f . This process continues 
until N c chromosomes have been selected from X". At the conclusion of the tournament selection 
process, some of the original chromosomes within Xy may not have been selected even once, while 
other chromosomes may have been selected multiple times. This selection algorithm is implemented by 
setting ISELECT = 2. 

Bin selection (version 1)— The first bin selection algorithm is different from greedy and tournament 
selection, as implemented here, in two general ways. First, the bin selection algorithm chooses its 
chromosomes from the accumulation file, not the active file. Second, the bin selection algorithm divides 
the distance along the Pareto front into N bin equal segments or “bins” using an arc-length computation 
and then selects an equal number of chromosomes from each segment. This ensures an equal 
distribution of selected points along the Pareto front. The parameter N bjn is a user-controlled parameter 
typically equal to 5 or 1 0. 

More precisely, this bin selection algorithm (version 1) is implemented as follows: 

(1) Initiate optimization using greedy selection and proceed until the total number of points on the Pareto 
front equals or exceeds N tot , which is a user controlled parameter that is typically several times the value 

N bin ■ 

(2) The arc-length along the Pareto front is computed, then divided into N bin equal arc-length segments or 
bins. Each existing Pareto front solution point is assigned to a bin according to its arc-length value. 


j = 1 N c 


it = \...,N c 


( 9 ) 
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Depending upon the distribution of points along the Pareto front, some bins may be heavily populated and 
others lightly populated. 

(3) The selection process randomly chooses a chromosome from the accumulation file keeping track of 
which bin it was chosen from. When the number of chromosomes chosen from a particular bin equals 
N avg defined by 


N avg N c / N bj n 

further selections from this bin are blocked. The process continues until N c chromosomes have been 
selected. For situations when N c is not evenly divisible by N bin , the value of N avg is incremented by one. 
Thus, certain bins may supply one more element to the next generation than other bins, but generally, the 
distribution will be reasonably well distributed across the entire Pareto front. This baseline variation of the 
bin selection algorithm (version 1) is implemented by setting ISELECT = -3. 

This selection algorithm does not guarantee that the Pareto front endpoints will be selected and carried 
forward to the modification process or passed through to the next generation. Since endpoint retention 
can be an attractive feature for selection, a variation of this selection algorithm that first selects the Pareto 
front endpoints and then proceeds with the standard selection algorithm is also available and is 
implemented when ISELECT = 3. 

Since the present bin selection algorithm requires the computation of arc-length along the Pareto front, it 
inherently requires optimization problems with two objectives, that is, the Pareto fronts must be a curve in 
two-dimensional space. Another bin selection option that works for two or more objectives is available and 
will be discussed next. 

Bin selection (version 2)— The second bin selection algorithm is similar to the first in philosophy and 
has virtually identical convergence properties, but differs in how the bins are constructed. It does not use 
an arc-length computation along the Pareto front. Instead, the design space is divided into a series of 
equally sized boxes or bins, each with N 0 -dimensions. This approach is applicable for any number of 
objectives (greater than one) and thus is more general than the first bin selection algorithm. 

More precisely, this bin selection algorithm (version 2) is implemented as follows: 

(1) Initiate the optimization using greedy selection and proceed until the total number of points on the 
Pareto front equals or exceeds N tot —defined the same as in bin selection algorithm version 1 . 

(2) The design space bins are constructed next by subdividing each objective dimension in the design 
space into M se0 equal segments— where M se0 is a user-controlled parameter. Thus, the design space is 

divided into a total of (M seg ) N ° bins. For example, when M seg =5 and N 0 =2, the design space is 
divided into 25 bins, each with the same rectangular shape. 

(3) Each existing element on the Pareto front, that is, each element in the accumulation file, is categorized 
as to which bin it occupies. Many of the design space bins will be empty, because they lie beyond the 
Pareto front or because they occupy an interior portion of the design space. Some of the bins that contain 
Pareto front elements may have many entries— some may have few. This, of course, depends on how the 
Pareto front is populated and how the Pareto front slices through the matrix of bins that have been 
created. The number of bins that have at least a single entry is counted. That resulting value is then 
stored in N bin . 

(4) The selection process randomly chooses a chromosome from the accumulation file keeping track of 
which bin it was chosen from. When the number of chromosomes chosen from a particular bin equals 
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N avg — defined the same as above— further selections from this bin are blocked. The process continues 
until N c chromosomes have been selected. As above, for situations when N c is not evenly divisible by 
N bjn , the value of N avg is incremented by one. Thus, certain bins may supply one more element to the 

next generation than other bins, but generally, the distribution will be reasonably well distributed across 
the entire Pareto front. This baseline variation of the bin selection algorithm (version 2) is implemented by 
setting ISELECT = -4. 

As with the bin selection algorithm (version 1), this selection algorithm does not guarantee that the Pareto 
front endpoints will be selected and carried forward to the modification process or passed through to the 
next generation. A variation of this algorithm that first selects the Pareto front endpoints and then 
proceeds with the baseline selection algorithm is also available and is implemented when ISELECT = 4. 

Modification Operators 

P Vector— After the new chromosomes have been selected and placed in the temporary holding array, 
G l , they must be modified using one of several modification operators to obtain the (n + 1) st generation of 
chromosomes, G n+1 . In the present implementation four modification operators are used— passthrough, 
random average crossover, perturbation mutation and mutation. How many chromosomes are modified 
with each operator is controlled by the P vector, which consists of four parameters— p B , p A , p P , p M . 
Each element of the P vector controls one modification operator. The value of each P vector element 
ranges from 0 to 1.0, and, for consistency, the sum of all four elements must always equal one. A P 
vector equal to 0.1, 0.3, 0.3, 0.3, for example, will cause the first 10 percent of the chromosomes to be 
modified using the passthrough operator, the next 30 percent to be modified using random average 
crossover, the next 30 percent to be modified using perturbation mutation and the last 30 percent to be 
modified using mutation. That is, p B =0.1, p A = 0.3, p P = 0.3, and p M = 0.3. 


The passthrough operator is always performed first. After passthrough is complete, the implementation 
order of the remaining operators is immaterial. Once all values of G n+1 have been established, the 
algorithm proceeds to fitness evaluation, ranking, and onto succeeding generations until the optimization 
is sufficiently converged. 

Passthrough— The simplest modification operator used in the present GA is “passthrough.” As the name 
implies, a certain number of chromosomes with the highest individual fitness values are simply “passed 
through” to the next generation from G f to G n+1 without modification. The number of chromosomes that 
are passed through to the next generation is controlled by the first parameter in the P vector, p B . The 

passthrough operator is always performed on the first p B N c chromosomes in G f . Care must be taken 
when choosing p B and N c so that p B N c >N 0 . If this is done and if a selection option that retains end 
points is used, the chromosomes with the highest individual fitness values will always be passed through 
to the next generation, thus guaranteeing that none of the individual maximum fitness values will ever 
decrease during GA iteration. 


Random average crossover— The next GA modification operator is called random average crossover 
and is implemented by first selecting two random chromosomes Xj-, and Xj 2 from G f . Next, the two 
selected chromosomes are combined on a gene-by-gene basis using the following formula: 


n + 1 
C 'J 


= ^( x /ji +x /, 72 ) i = X2,...,N e 


( 10 ) 


where x"^ 1 corresponds to the / th gene in the y th chromosome associated with G n+1 and x,- ;1 andx/ ;2 
correspond to the / th genes from the randomly chosen chromosomes X^ and Xy 2 . The number of 
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chromosomes modified using the random average crossover operator is determined by the parameter 
p A — the second element in the P vector. 

Perturbation mutation— The next GA modification operator is called perturbation mutation and is 
implemented by first selecting a random chromosome Xy from G f . Next, a probability test is performed 
for each gene xj j in the selected chromosome involving a call to the random number generator R(0,1). If 
the returned random number is less than p 1 the gene is modified using 

= xjj + (x max . -* min ,)[fl(0 I 1)-0.5]/3 (11) 

where (5 is a user-specified tolerance that controls the size of the perturbation mutation, and p 1 is a user- 
specified control parameter that statistically controls the number of genes that are modified. For sensible 
results the values of /3 and p, must be between 0 and 1 .0. 

Because this operator can cause the value of a particular gene to exceed one of its constraints (x max or 
x min ), checks are required to prevent this. The number of chromosomes modified using the perturbation 
mutation operator is determined by the parameter p P —the third element in the P vector. 

Mutation— The last GA modification operator used in the present study is called mutation and is 
implemented similarly to the perturbation mutation operator. First, a random chromosome X' is chosen 
from G l . Next, a probability test is performed for each gene xjj in the selected chromosome involving a 
call to the random number generator R(0,1). If the returned random number is less than p 2 the gene is 
given a completely different value using 

x if=(x max, -*min,)fl(0.1)+X min , (12) 

The parameter p 2 is a user-specified control parameter that statistically controls the number of genes that 
are modified. For sensible results p 2 must be between 0.0 and 1.0. The number of chromosomes 
modified using the mutation operator is determined by the parameter p M —the fourth element in the P 
vector. 

Gene-space transformation — An option for accelerating GA convergence for multi-objective 
optimization problems involving a gene-space transformation is described next. This option, introduced in 
Holst and Pulliam, 15 worked amazingly well for model problems with simple non-convoluted Pareto fronts, 
but performed poorly for model problems in which the Pareto front was convoluted. Thus, it is of interest 
to see how the gene-space transformation procedure will perform on a real-world problem in which the 
shape of the Pareto front— in gene space — is unknown. 

The transformation algorithm is outlined as follows: 

(1) Transform G l using the gene-space transformation procedure. 

(2) Perform modifications on the transformed chromosomes according to the user-input P vector 
values. 

(3) Using the inverse of the original transformation the modified chromosomes are transformed back 
to obtain G n+1 . 

This gene-space transformation procedure, which can be viewed as a gene-space rotation of coordinates, 
causes a linear coupling between each of the genes, and thus, affects how they are changed in the 
modification process. In some cases the gene-space transformation procedure can significantly improve 
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GA convergence. The transformation procedure, which theoretically works for any number of objectives, 
will now be presented for the two-objective case, that is, for N 0 = 2. 

The idea behind the transformation procedure is to perform a rotation of coordinates in gene space using 
an angle-preserving, length-preserving orthogonal transformation. For this purpose a simplified version of 
the Gram-Schmidt orthogonalization is used. 45 The current set of n th generation chromosomes (those that 
have been newly selected and placed in the temporary holding array G f ) can be written as 

X 1,2 ■■■ x \n c 

i 

X N S ,N C 



where each column is one of the chromosomes in the active file holding array. It is this matrix that is to be 
transformed using 

UGi = G n (13) 

where U is a unitary matrix of rank N G that needs to be constructed and G n is the resulting transformed 
matrix. To construct U a set of basis vectors R needs to be established. The elements of R are defined 
by 


/ /,y=X 2 -X 1 for j = 1 


(14a) 


and 


f 1 if i = j 

r, : = for / > 2 

IJ [0 if i + j 1 


The simple choices made for r, y when j > 2 serve to simplify the transformation matrix construction 
without sacrificing overall generality. 

As can be seen from Eq. (14a), the first coordinate direction in the new transformed coordinate system is 
chosen to be parallel to X 2 -X v As was mentioned earlier, the chromosome with the best fitness for the 
first objective is always placed into X^ and the chromosome with the best fitness for the second objective 
is always placed into X 2 . This convention is crucial for the success of the transformation algorithm, as it 
causes the first transformation coordinate to be aligned with the Pareto front endpoints. Keep in mind that 
it is imperative that the first element of the P vector— the element that controls passthrough — is large 
enough so that both X., and X 2 are always passed through to the next generation without modification. It 
is also imperative for the selection algorithm to have an endpoint retention option. 

The unitary transformation matrix U is constructed column by column using 

first column 

u it i=— (15a) 

n. i 
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second column 


Vi-n,2-u 2 ,M,i ■ 



n th column 


Vi = n, n 


/7 — 1 

- S U n,j u i,j > 


y-i 



(15b) 


(15c) 


Once G n is obtained the modification operators are applied the same as without transformation to obtain 
G n+1 . The final G n+1 values are then obtained using an inverse transformation indicated by 

U r G n+1 = G n+1 (16) 

Because U is a unitary matrix, its inverse is simply its transpose, and thus, it is easily constructed. 

The gene space transformation option is controlled using the ITRAN control parameter. If ITRAN = 0, no 
gene space transformation in implemented. If ITRAN = 1, the gene space transformation algorithm is 
activated. 


Geometry Parameterization 

In aerodynamic shape optimization, the geometric parameterization is an important step that effectively 
connects the aerodynamic analysis routine to the optimization routine. An analytic parameterization 
suitable for wing and wing-fuselage configurations, typical of the transonic flow regime, is now described. 
Most of the parameters have common-sense definitions in which the name itself provides the definition, 
for example, the wing leading edge radius of curvature, the wing leading edge sweep, or the fuselage 
length. 

A Cartesian coordinate system (x,y,z ) is used for all configurations. The coordinate system origin is at 
the wing root leading edge for isolated wing configurations and at the fuselage nose for wing-fuselage 
configurations. A plane of symmetry (y = 0) is inherently assumed for all configurations. The x 
coordinate is aligned with the freestream direction and is positive downstream, y is normal to the 
symmetry plane, positive to the pilot’s right, and z is vertical, positive upward. All length parameters are 
nondimensionalized using the wing-root chord, that is, the wing chord length where it intersects the plane 
of symmetry. 

Isolated-wing configurations — Isolated-wing geometries are parameterized using N w airfoil defining 
sections or stations, each using the geometric parameterization of Sobieczky. 46 The first defining station is 
always at the wing root and the last is always at the wing tip. Each defining airfoil section is characterized 
by ten parameters. These parameters along with brief definitions are listed in Table 1. A graphical 
description of these parameters is presented in Fig . 1 . All angles are measured in degrees. Once these 
parameters are specified the airfoil coordinates (z as a function of x) are determined using a polynomial 
of the form 


z = 


2 > n x "- 1/2 


n = 1 


(17) 


where the coefficients a n are computed from the ten airfoil defining parameters. Two applications of Eq. 
17 are required for a complete airfoil specification— one for the upper surface coordinates and one for the 
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lower surface coordinates. Six simultaneous linear equations involving rle, xup, zup, zxxup, 
ate + /> te/2, and zte are solved for the upper surface coordinates, and six simultaneous linear equations 
involving rle, xlo, zlo, zxxlo, ate- /3te/2 , and zte are solved for the lower surface coordinates. Because 
rle and zte are used for both surfaces, slope continuity at the leading edge and zero thickness at the 
trailing edge are always maintained. 



AXIAL COORDINATE - x/c 


Fig. 1 Airfoil parameterization used for each wing defining station (taken from Sobieczky 46 ). 
Angles are measured in degrees. All lengths are nondimensionalized by wing-root chord. 

Once the airfoil coordinates are constructed for each defining station, wing coordinates are constructed 
using linear lofting and then modified according to a set of planform parameters. Seven planform 
parameters are utilized at each airfoil defining station, producing a total of 7xN w planform parameters. 
These parameters are also listed in Table 1 along with brief definitions. 

Certain geometric parameters have predetermined values or are not formerly used in the wing definition 
process. For example, the root chord length is always unity, and the spanwise location for the root airfoil 
defining station is always zero. The tip airfoil defining station values for the dihedral angle and the leading 
edge sweep angle have no meaning, but are nevertheless included in the list for each wing 
parameterization. 


Table 1. Definitions for the airfoil section and wing planform parameters. A value for each 
parameter is required at each of the N w wing defining stations, j = N w . 


Symbol 

Definition 

Symbol 

Definition ! 

Airfoil section parameters 

Planform parameters 

rle, 

Leading edge radius of curvature 

c i 

Airfoil chord length 

xup. 

x-location of upper surface max 
thickness 

D j 

Dihedral angle (deg) 

zup, 

Upper surface max thickness 

A ; 

Leading edge sweep (deg) 

zxxup, 

Upper surface curvature at x =xup, 


Airfoil section thickness multiplier 

xlo, 

x-location of lower surface max 
thickness 

d i 

Wing twist (deg) 

zlo, 

Lower surface max thickness 

x e,j 

Center of twist 

zxxlo, 

Lower surface curvature at x =xlo , 

Vi 

Spanwise defining station location 

zte,- 

Position of trailing edge above z = 0 



ate, 

Angle between trailing edge bisector 
and z = 0 plane (deg) 



e, 

Trailing edge included angle (deg) 




15 






Wing-fuselage configurations— Wings for wing-fuselage configurations are parameterized using the 
procedure as described above. Thus, only the fuselage portion of the wing-fuselage parameterization will 
be discussed in this section. The first step is to define a base fuselage, which is analytically constructed in 
three sections. The nose section is an ellipsoid of revolution. The main body section is a right circular 
cylinder that smoothly transitions into the nose section. And the boattail section is a sine curve of 
revolution smoothly connecting the main portion of the fuselage with a downstream sting. The sting — not 
formally part of the fuselage — is a small-radius right circular cylinder that extends to the downstream 
outflow boundary. A total of eight parameters— defined in Table 2— are required to specify this base 
fuselage. Note: The y-location of the wing-root leading edge is always assumed to be on the symmetry 
plane. Thus, y R is inherently assumed to be zero and is only included for completeness. 

Table 2. Base fuselage parameter definitions. All lengths are nondimensionalized by airfoil root 
chord. 


Symbol 

Definition 

X R 

x-location of wing root leading edge 

y r 

y-location of wing root leading edge 

z R 

z-location of wing root leading edge 

X B 

Body length 

X N 

Length of ellipsoid nose 

r B 

Radius of cylindrical portion of 
fuselage 

X T 

Length of fuselage boattail 

r s 

Radius of fuselage sting 


Modification of the base fuselage shape is achieved using a series of Hicks-Henne bump functions. 47 A 
total of N x bump functions are distributed axially with equal spacing along the x-direction for each of N T 
circumferential stations also distributed with equal spacing in the circumferential direction. This allows a 
total of N x xN T decision variables that are used to produce incremental radius perturbations Ar's to the 
base fuselage. Any position on the base fuselage, except the nose x = 0 or tail x =x B , can be changed 
using this approach. A sketch showing the definitions of these parameters, including one possible 
arrangement of Ar's is shown in Fig. 2. 



Fig. 2 Sketch showing decision variable definitions for the fuselage parameterization. For this 
arrangement there are three axial stations ( N x =3) and four circumferential stations (N t = 4). 

Masking array— Associated with each decision variable, that is, each gene, is an integer value— either 
one or zero— stored in a masking array called mask . The masking array is established at the beginning 
of the optimization computation and tells the GA which parameters to modify in the optimization 
process— using the crossover and mutation operators— and which to leave unmodified. For example, the 
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mask values for c-,, y v A Nw , D w ^and y R are always zero. These parameters, as described above, are 

carried along with each chromosome but are never used by the GA process. The masking array can also 
be used to limit the size of the design space. Only those genes with mask = 1 are utilized in the search for 
optimal problem objectives. 


Computed Results 

Area Error Norm for Two-Objective Problems 

In order to evaluate the accuracy and efficiency of an optimization algorithm it is important to have a 
proper error norm to assess the level of convergence. For single objective problems a suitable error norm 
is the simple arithmetic difference between the current “best fitness” and the exact answer. This, of 
course, works well for model problems in which the exact answer is known. For applications in which the 
exact answer in not known a priori, it can still be employed as an a posteriori error computation. 

For multi-objective optimization, a workable error norm is more elusive since a range of solution values — 
the so-called Pareto front— is being sought. The topic of how to define easy-to-implement and accurate 
error norms for comparing one Pareto front approximation to another for the same problem is discussed 
at length in Zitzler, et al. 11 and Knowles and Come. 18 A suitable norm encompassing the optimal values of 
each of the individual objectives is one possibility. Flowever, such a norm would not take into account the 
trade-off regions of the Pareto front. Another possibility is the attainment surface approach of Fonseca 
and Fleming. 48 In this approach a set of equally spaced sampling lines that intersect the full breadth of the 
Pareto front are used. Statistical measures of goodness can then be developed based upon how many 
intersection points one Pareto front has that are superior to another Pareto front. 

In the present study the area between the current Pareto front and the exact Pareto front— a numerically 
computed approximation of the exact Pareto front— is used as the error norm to determine level of 
convergence. As the size of this quantity goes to zero, the current solution approaches the final solution. 

The area is computed numerically by dividing the region between the “exact” Pareto front and the current 
approximation into a series of triangles as shown in Fig. 3. The area error norm is the sum of each of the 
individual triangular areas. Thus, an approximation to the exact Pareto front that fails to match the final 
solution in any location will produce a non-zero value for the area error norm. 



Fig. 3 Area error norm computational process for a two-objective Pareto front. 
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Isolated wing results 


The first computed result involves a two-objective optimization for an isolated wing flying in the transonic 
speed regime— =0.84. The two objectives — lift-to-drag ratio at fixed lift and the inverse of the 
nondimensional structural mass— are simultaneously maximized by the previously presented MOGA 
optimization procedure. The aerodynamic lift-to-drag ratio is evaluated for each candidate design, that is, 
for each chromosome, using the TOPS full potential solver 49 by iterating on angle of attack to obtain the 
lift-to-drag ratio at the specified value of lift— C L = 0.45. The tolerance on the lift iteration is ±1%, but in 
most cases the error is much less. In effect, this is drag minimization at fixed lift. 

The structural mass computation utilizes an equivalent box beam model similar to that used in Oyama. 42 
Given a set of loads from the aerodynamics routine, the structures model computes the minimum-mass 
box beam that will support that load with a factor of safety of two. Shear and bending are included in this 
model but not torsion. 

The design space parameterization for isolated wings with two defining stations, N w = 2, as described 
above, consists of 34 parameters. Airfoil defining station one is at the wing root, and station two is at the 
wing tip. Maximum and minimum constraint values for each of the geometric parameters, as well as 
masking array values, are given in Table. 3. Except for wing twist, all planform geometric parameters 
have a masking array value of zero, that is, they are not modified in the optimization process. Thus, a 
total of 22 genes are active within each chromosome for this optimization problem, that is, N G =22. 

Table 3. Maximum and minimum gene constraints along with masking array values for the 
isolated win g optimization problem. 


Parameter 

max 

min 

mask 

Parameter 

max 

min 

mask 

Defining airfoil section 1 

Defining airfoil section 2 

rle! 

0.02 

0.008 

1 

rle 2 

0.02 

0.008 

1 

xup! 

0.45 

0.35 

1 

xup 2 

0.45 

0.35 

1 

zupi 

0.07 

0.06 

1 

zup 2 

0.07 

0.06 

1 

zxxup 1 

-0.2 

-0.8 

1 

zxxup 2 

-0.2 

-0.8 

1 

XlO-, 

0.45 

0.35 

1 

xlo 2 

0.45 

0.35 

1 

ZlO-, 

-0.06 

-0.07 

1 

zlo 2 

-0.06 

-0.07 

1 

zxxlo-, 

0.8 

0.2 

1 

zxxlo 2 

0.8 

0.2 

1 

zte.. 

0.01 

-0.01 

1 

N 
1 1 
CD 

N> 

0.01 

-0.01 

1 

ate! 

8.0 

0.0 

1 

ate 2 

8.0 

0.0 

1 

/3tei 

15.0 

6.0 

1 

P te 2 

15.0 

6.0 

1 

Planform defining station 1 

Planform defining station 2 

Cl 

1.0 

1.0 

0 

c 2 

0.5 

0.5 

0 

Di 

0.0 

0.0 

0 

d 2 

0.0 

0.0 

0 

A, 

37.0 

37.0 

0 

a 2 

37.0 

37.0 

0 

fi 

1.0 

1.0 

0 

^2 

1.0 

1.0 

0 


3.0 

0.0 

1 

0 2 

0.0 

-3.0 

1 

x e,i 

0.5 

0.5 

0 

x e,2 

0.5 

0.5 

0 

y i 

0.0 

0.0 

0 

y 2 

2.4 

2.4 

0 


Results from this multi-discipline optimization problem are presented in Figs. 4-9. For all computations in 
this series of results A/ c =32, /3 = 0.1,p 1 = 0.2, p 2 = 0.2 and P = (0.04,0.32,0.32,0.32). This P vector 
resulted in two passthrough chromosomes and ten chromosomes in each of the three remaining 
modification categories — random average crossover, perturbation mutation, and mutation. 
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Pareto front comparisons are presented in Figs. 4. Note that each comparison contains five curves that 
correspond to different ISEED initialization values used in the random number generator. There is also a 
“master” result plotted in each part of Fig. 4 that is a concatenation of all Pareto front data produced for 
this case— a total of twenty curves— with only the non-dominated points retained. 



a) ISELECT = 1, ITRAN = 0 



b) ISELECT = 1, ITRAN = 1 



LIFT TO DRAG RATIO - L/D 


c) ISELECT = 3, ITRAN = 0 



d) ISELECT = 3, ITRAN = 1 


Fig. 4 Pareto front variation for transonic wing optimization in which the lift-to-drag ratio at fixed 
lift and the inverse of the structural mass are simultaneously maximized. Five solutions utilizing 
different seed values are plotted for each GA variation after 500 generations: = 0.84, 

C L = 0.45 , N c =32 N G =22, /3 = 0.1, p^= 0.2, p 2 =0.2, and P = (0.04, 0.32, 0.32, 0.32). 
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There are several interesting features to note in this series of Pareto front comparisons: 

First, the variation in level of convergence across the various ISEED values is quite striking, especially for 
the cases involving greedy selection without gene space transformation, ISELECT = 1, ITRAN = 0 [Fig. 
4(a)]. For these cases, the optimal lift-to-drag ratio point on the Pareto front converges consistently in the 
number of generations allotted, regardless of the ISEED value, but the minimum structural mass point 
does not converge consistently. 

This situation is improved, somewhat, for the ISELECT = 1, ITRAN = 1 case [Fig. 4(b)] where the best 
structural mass endpoint values are obtained. However, this latter set of curves suffers dramatically in the 
middle portion of the Pareto front where the two objectives form trade-offs with each other. 

The best results for this multi-discipline isolated-wing optimization are realized when bin selection is 
utilized, ISELECT = 3 [Figs. 4(c) and 4(d)]. The Pareto fronts for both sets of curves that utilize bin 
selection more closely approximate the master curve than their counterparts using greedy selection. The 
set of curves in Fig. 4(c)— the ISELECT = 3, ITRAN = 0 case— are overall the best in that they 
approximate the master curve most closely and have the smallest amount of statistical variation. 
Nevertheless, the minimum structural mass point for this set of computations still experiences a moderate 
amount of inconsistency in convergence. 

Convergence history comparisons for each of the Pareto front curves presented in Figs. 4 are presented 
in Fig. 5. Computation of the Pareto front error was achieved using the area error norm described above. 
The master Pareto front displayed in Figs. 4 was used as an approximation to the “exact” solution. In 
addition, for each set of convergence history curves displayed in Figs. 5, there is an average convergence 
history curve plotted (the dashed curve), which is computed by arithmetically averaging each point in the 
convergence history data file. 

Note that occasionally the convergence history error increases with generation. This usually occurs early 
in the convergence process when the Pareto front is still crudely defined with a small number of discrete 
points and is explained in the series of sketches shown in Figs. 6. Figure 6(a) shows a hypothetical two- 
objective Pareto front, both the exact curve and a crude approximation to it that might exist after the 
search has proceeded for n generations. Figures 6(b) and 6(c) show two different scenarios for what the 
Pareto front might look like after n + 1 generations. 

In both scenarios a single new point has been added which neither dominates nor is dominated by any of 
the previously existing points on the Pareto front. In scenario 1 the new point’s placement is less 
advantageous than that of scenario 2. Nevertheless, both new points are possible outcomes of the 
genetic algorithm search process. Despite the fact that scenario 1 represents an “enrichment” of the 
Pareto front, the area error norm actually increases in going from the n th generation to the n + 1 st 
generation. This is why the area error norm occasionally jumps during solution evolution. 

There is a moderately large amount of variance in each of the convergence history curves displayed in 
Fig. 5. This emphasizes the point that genetic algorithms are stochastically based. Comparison of results, 
especially convergence efficiency results, must be made with the proper amount of statistical averaging in 
place. Results from each MOGA variation demonstrate this conclusion. 
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a) ISELECT = 1, ITRAN = 0 


c) ISELECT = 3, ITRAN = 0 




b) ISELECT = 1, ITRAN = 1 


d) ISELECT = 3, ITRAN = 1 


Fig. 5 Convergence efficiency comparisons for transonic wing optimization in which the lift-to- 
drag ratio at fixed lift and the inverse of the structural mass are simultaneously maximized. Five 
solutions utilizing different seed values are plotted for each GA variation: 0.84, C L = 0.45, 

N c =32 N G =22, = 0.1, p^= 0.2, p 2 = 0.2 and P = (0.04,0.32,0.32,0.32). 


The averaged convergence history efficiencies for each of the four MOGA variations studied herein are 
presented on the same set of axes in Fig. 7. Note that the two variations that utilize bin selection are both 
faster than their greedy selection counterparts, achieving more than an order of magnitude reduction in 
the area error norm after 500 generations. The gene-space transformation variation produced mixed 
results that were slightly faster when used with greedy selection and slightly slower when used with bin 
selection. 
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OBJECTIVE 1 



OBJECTIVE 1 



OBJECTIVE 1 


a) Baseline, Pareto front after 
n generations. 


b) Pareto front after n+1 
generations (scenario 1). 


c) Pareto front after n+1 
generations (scenario 2). 


Fig. 6 Sketch of Pareto front convergence history process showing two different advancement 
scenarios. 



GENERATION NUMBER 

Fig. 7 Average convergence history curves for GA variations plotted in Figs. 4 and 5: M x =0.84, 
C L = 0.45 , N c =32 N g =22, 0-0.1, p, - 0.2, p 2 = 0.2 and P = (0.04,0.32,0.32,0.32). 

Computed results for the two-objective isolated-wing optimization problem taken from the endpoints of a 
suitable ISELECT = 3, ITRAN =0 Pareto front are displayed in Figs. 8-10. Figure 8(a) shows the upper 
surface Mach number contours for the minimum drag point (at fixed lift), and Fig. 8(b) shows the same 
view for the minimum structural mass point. Note the extreme suction pressure (high Mach number) 
upstream of the shock wave in Fig. 8(b). This is caused by a large amount of flow expansion on the 
forward part of the wing and results in a much stronger shock wave than the solution depicted in Fig. 8(a). 

This situation can be viewed in a more quantitative way by looking at cross-sectional plots of the pressure 
coefficient, which are displayed in Figs. 9 for three spanwise stations — y / b = 0.21, 0.51, 0.78. In each plot 
the pressure coefficient for the two Pareto front endpoints — minimum drag and minimum structural 
mass— are displayed. This series of figures shows that the minimum drag case (solid curves) has a 
greatly reduced shock strength relative to the minimum structural mass case (dashed curves). 

Both solutions have the same total lift, but the spanwise distribution of lift is different between the two 
cases. This can be seen from Figs. 9 by carefully examining the areas circumscribed by each cross- 
sectional pressure coefficient plot. The minimum structural mass solution (dashed curves) contains more 
area inboard and less outboard relative to the minimum drag point. This tends to move the center of 
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pressure inboard for the minimum structural mass solution, and thus, reduces the moment arm through 
which the wing lift acts. This, in turn, reduces the wing root bending moment and allows for a lighter 
structure. 



a) Min drag b) Min structural mass 


Fig. 8 Upper surface Mach number contours for transonic wing optimization in which drag at fixed 
lift and structural mass are simultaneously minimized: ISELECT = 3, ITRAN = 0, M x = 0.84, 
C L = 0.45 , N c = 32 N G =22, 0-0.1, ^=0.2, p 2 = 0.2 and P = (0.04,0.32,0.32,0.32). 





Fig. 9 Pressure coefficient distributions at selected spanwise stations for the wing optimization 
problem presented in Fig. 8 showing differences between the two Pareto front endpoints: 

ISELECT = 3, ITRAN = 0, 0.84, C L =0.45, N c =32 N G =22, 0 = 0.1, ^=0.2, p 2 = 0.2 and 

P = (0.04, 0.32, 0.32, 0.32). 

Another important attribute in this multi-discipline optimization problem is wing thickness— especially at 
the wing root. Selected airfoil cross-sectional profiles are presented in Fig. 10 at y/b = 0.0,0.48,0.98 
(wing root, mid-span, wing tip). Profiles for both Pareto front endpoints — minimum drag (solid curves) and 
minimum structural mass (dashed curves)— are presented. Note that the vertical axis had been greatly 
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expanded to facilitate viewing of the results, and that each profile is displayed in rigged position with wing 
twist included. For this set of computations the minimum-drag root-airfoil thickness is the minimum 
allowed by the geometric constraints, and the minimum-structural-mass root-airfoil thickness is the 
maximum allowed. This is as expected. A reduced wing-root thickness helps reduce drag by reducing 
shock strength and thus reduces wave drag. An increased wing-root thickness helps to reduce structural 
mass by more efficiently supporting the wing-root bending moment with a box beam that has increased 
depth. 



NORMALIZED AXIAL DISTANCE - x/c 

Fig. 10 Selected airfoil profiles in rigged positions— at y/b=0.0, 0.48 and 0.98— for the transonic 
wing optimization problem of Fig. 8 showing differences between the two Pareto front endpoints 
(note the expanded z/c scale): ISELECT = 3, ITRAN = 0, M ao = 0.84, C L = 0A5, N c =32, N G =22, 
P = 0.1, Pi = 0.2 , p 2 = 0.2, and P = (0.04,0.32,0.32,0.32). 

Wing-fuselage results 

The next set of computed results involves a two-objective optimization for a wing-fuselage configuration 
flying in the transonic flow regime— =0.84. The two objectives— lift-to-drag ratio at fixed lift and the 
configuration’s volume — are simultaneously maximized by the previously presented MOGA optimization 
procedure. As before, the aerodynamic lift-to-drag ratio is evaluated for each candidate design, that is, for 
each chromosome, using the TOPS full potential solver 49 by iterating on angle of attack to obtain the lift- 
to-drag ratio at the specified value of lift— C L = 0.45 . 

The wing-fuselage volume computation utilizes a straightforward algorithm that divides the fuselage into a 
series of cross-sections, computes the area of each cross-section using a simple triangularization, and 
then computes the volume by integrating along the fuselage using trapezoid integration. The wing volume 
is computed using a similar algorithm that is deactivated for those wing cross-sections that lie inside the 
fuselage. 

There are several reasons lift-to-drag ratio (at fixed lift) and wing-fuselage volume have been chosen as 
the two objective functions for the present optimization problem. First, there is a clear-cut conflict between 
these two objectives. When the volume is maximized, the lift-to-drag ratio suffers and vice versa. Thus, 
we have a classical engineering trade-off problem. Flow the MOGA performs on such a problem with real- 
world objectives is the key aspect of this study, not the development of new aeronautical engineering 
concepts. 


24 




In addition, the vehicle volume was chosen as an objective because it is an easy quantity to understand. 
The direct computation of maximum volume is easy to perform, thus providing a simple check for the 
results obtained by the GA— at least for one endpoint on the Pareto front. This problem is interesting in 
that it is complicated by transonic-flow shock-wave-induced nonlinearities and by a large number of 
decision variables that have large sensitivity variations. 

The present wing-fuselage parameterization with two wing defining stations, N w = 2, as described above, 
consists of 66 parameters. Airfoil defining station one is at the wing root, and station two is at the wing tip. 
There are 24 A r values used to define the fuselage— six equally-spaced axial stations, each with four 
equally-spaced circumferential locations. The axial positions are located at 
x/c = 0.71, 1 .43, 2.14, 2.86, 3.57, 4.29 where the fuselage length is fixed, x B =5.0. The circumferential 
positions are located at cp = - 90°, -30°, 30° and 90° where cp = - 90° corresponds to the fuselage keel 
line, and cp =90° corresponds to the crown line. 

Maximum and minimum constraint values for each of the geometric parameters, as well as masking array 
values, are given in Table. 4. As can be seen, 43 of the 66 genes that make up each chromosome have 
masking array values of one, and 23 have values of zero, that is, only 43 genes are actually modified 
during the optimization process — N a =43. Note also that the max-min values associated with each of the 
A r values are not symmetric about zero. Along the keel and crown lines, for example, any A r value 
averaged between the max and min constraints, will be positive, and along the other meridinal stations 
(along the fuselage side), the same averaging will produce negative values. The max and min A r 
constraints were chosen in this way to keep the vertical extent of the fuselage from becoming too small. A 
difficulty arises in the wing-fuselage line of intersection computation if the intersection line contacts the 
symmetry plane. This results in fuselage designs that have depth-to-width ratios that generally exceed 
unity— sometimes by a substantial amount. 

The first wing-fuselage optimization results are presented in Figs. 11-15. For these results the freestream 
Mach number is 0.84 and the lift coefficient is held fixed at 0.45. The tolerance on the lift iteration is ±1%, 
but in most cases the error is much less. The second binning selection algorithm with endpoint retention 
is used, ISELECT = 4, and the gene-space transformation procedure is turned off, ITRAN = 0. Both 
mutation probabilities are fixed at 0.2 and P = (0.04, 0.32, 0.32, 0.32) . 

This optimization uses 34 chromosomes for each generation. The first element of the P vector is set such 
that two passthrough chromosomes are chosen during each generation— always the Pareto front 
endpoints— the current minimum drag and the maximum volume chromosomes. Thus, 32 modified 
chromosomes require new function evaluations during each generation. These function evaluation 
computations are efficiently performed simultaneously on 32 processors of a parallel computer. The other 
three elements of the P vector are set to evenly divide the remaining chromosomes between the three 
remaining modification operators — random average crossover, perturbation mutation, and mutation (at 
least as close to even as possible). 

Development of the Pareto front as a function of generation number (n) is displayed in Fig. 11 for two 
cases— a baseline case and a second case in which the mutation process was aggressively biased 
toward the design space boundary. For the baseline case each mutated gene value was chosen to lie 
between that particular gene’s maximum and minimum constraints using a random process, as previously 
described. 

For the second case — hereafter called the aggressive mutation case— the perturbation and global 
mutation operators were modified. Each new gene value was randomly chosen from a much larger 
range — an order of magnitude larger than the range established by the user-specified constraints. If the 
resulting gene value fell between the original constraints for that gene, it was used without modification. If 
it was greater than the maximum constraint or less than the minimum constraint, it was set to the nearest 
constraint limit. All other aspects of the MOGA approach as outlined above were unaltered. 
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Table 4. Maximum and minimum gene constraints along with masking array values for the wing- 


f uselage optimization problem. 


Parameter 

max 

min 

mask 

Parameter 

max 

min 

mask 

Parameter 

max 

min 

mask 

Defining airfoil section 1 

Defining airfoil section 2 

Fuselage perturbations 

rle! 

0.014 

0.014 

0 

rle 2 

0.014 

0.014 

0 

A/v, 

0.05 

-0.02 

1 

xup! 

0.45 

0.35 

i 

xup 2 

0.45 

0.35 

i 

A/2,1 

0.05 

-0.02 

1 

zup. 

0.08 

0.05 

i 

zup 2 

0.08 

0.05 

i 

A/3,1 

0.05 

-0.02 

1 

ZXXUP! 

-0.2 

-0.8 

i 

zxxup 2 

-0.2 

-0.8 

i 

A/4,1 

0.05 

-0.02 

1 

XlO! 

0.45 

0.35 

i 

xlo 2 

0.45 

0.35 

i 

A/5,1 

0.05 

-0.02 

1 

Z\Oi 

-0.05 

-0.08 

i 

zlo 2 

-0.05 

-0.08 

i 

A/6,1 

0.05 

-0.02 

1 

zxxlo 1 

0.8 

0.2 

i 

zxxlo 2 

0.8 

0.2 

i 

A/), 2 

0.02 

-0.05 

1 

zte.. 

0.0 

-0.0 

0 

zte 2 

0.0 

-0.0 

0 

Ar 22 

0.02 

-0.05 

1 

ate! 

8.0 

0.0 

1 

ate 2 

8.0 

0.0 

1 

A/3,2 

0.02 

-0.05 

1 

P te! 

10.0 

10.0 

0 

/3te 2 

10.0 

10.0 

0 

A/4 2 

0.02 

-0.05 

1 

Planform defining station 1 

Planform defining station 2 

A/5,2 

0.02 

-0.05 

1 

Cl 

1.0 

1.0 

0 

c 2 

0.5 

0.5 

0 

Ar 6 , 2 

0.02 

-0.05 

1 

Di 

0.0 

0.0 

0 

d 2 

0.0 

0.0 

0 

A/), 3 

0.02 

-0.05 

1 

A, 

37.0 

37.0 

0 

a 2 

37.0 

37.0 

0 

CO 

< 

0.02 

-0.05 

1 


1.0 

1.0 

0 

^2 

1.0 

1.0 

0 

A/3,3 

0.02 

-0.05 

1 

o, 

3.0 

0.0 

1 

g 2 

0.0 

-3.0 

1 

> 

CO 

0.02 

-0.05 

1 

X B, 1 

0.5 

0.5 

0 

X B, 2 

0.5 

0.5 

0 

CO 

V? 

< 

0.02 

-0.05 

1 

y i 

0.0 

0.0 

0 

y 2 

2.4 

ea 

0 

CO 

< 

0.02 

-0.05 

1 

Base fuselage parameters 

A/1,4 

0.05 

-0.02 

1 

X R 

1.5 

1.5 

0 

X N 

0.5 

1.5 

1 

A/2,4 

0.05 

-0.02 

1 

y r 

0.0 

0.0 

0 

r B 

0.35 

0.35 

0 

A/3,4 

0.05 

-0.02 

1 

z R 

-0.15 

0.15 

1 

X T 

0.5 

1.5 

1 

A/4,4 

0.05 

-0.02 

1 

X B 

5.0 

5.0 

0 

r s 

0.1 

0.1 

0 

A/5,4 

0.05 

-0.02 

1 









A/6,4 

0.05 

-0.02 

1 


Development of the aggressive mutation approach was motivated by realizing that the maximum volume 
endpoint for the present problem lies on the gene space boundary. In particular, 38 of the 43 modifiable 
genes are maximized or minimized by the MOGA in achieving the maximum volume solution. In addition, 
the minimum drag point (at fixed lift) is achieved with 27 maximum or minimum gene values. Thus, a 
mutation process that biases toward the design space boundary seems to be logical and, for the present 
case, as seen by comparing Figs. 1 1 (a) and 1 1 (b), produces a significant acceleration in convergence. 

For the aggressive mutation case, 96% of the normalized maximum volume— which occurs at a value of 
1.416 — is achieved after 50 generations. The baseline MOGA solution, after 50 generations, has 
achieved only 85% of the maximum volume. The minimum drag point is not significantly affected by the 
mutation strategy, converging at about the same rate for both approaches. The convergence 
characteristics, displayed in Fig. 11 for one ISEED value, were generally observed for every ISEED value 
that was tested for the present wing-fuselage problem. 

For multi-objective optimization problems in which the Pareto front lies inside the gene-space interior, the 
aggressive mutation approach may not be advantageous. More research is needed to answer this 
question. All additional cases presented herein utilize the aggressive mutation approach. 
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a) Baseline 


b) Aggressive gene mutation 


Figure 11: Pareto front convergence for wing-body optimization involving drag minimization and 
volume maximization: C L = 0.45, M^ = 0.84, ISELECT = 4, ITRAN = 0, N 0 =2, A/ c =34, N G = 43. 


Several physical characteristics associated with the Fig. 11(b) design space are presented in Fig. 12. The 
Pareto front, identical to the n = 500 curve from Fig. 11(b), is duplicated as a solid line. Each dot lying 
below and to the left of the Pareto front represents one function evaluation (one call to the TOPS analysis 
code). The three positions called out along the Pareto front— A, B, and C — indicate the design-space 
locations for the three pressure distributions plotted on the right. They correspond to solutions at A) 
maximum volume, B) intermediate, and C) minimum drag (at fixed lift). 



Figure 12: Design space for wing-body two-point optimization problem: C L = 0A5, M CO =0.84, 
ISELECT = 4, ITRAN = 0, N 0 =2, N c =34, N G =43. Wing sectional pressures at three positions 
along Pareto front: A) max volume, B) intermediate, C) min drag. 
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For each solution two airfoil pressure distributions are presented— the first inboard near the wing-fuselage 
juncture ( y/b = 0.2) and the second near the mid-semi-span ( y/b = 0.6). The most obvious feature from 
these results is the change in shock strength as the Pareto front is traversed. The maximum volume 
solution (A) contains strong shocks on both the upper and lower wing surfaces, while the minimum drag 
solution (C) contains mild shocks on only the wing’s upper surface. 

Upper surface Mach number contours for each of the three solutions presented in Fig. 12 are displayed in 
Figs. 13. This series of figures shows the shock wave pattern variation on the upper surface of the wing- 
fuselage configuration as the solution transitions from the maximum volume Pareto front endpoint to the 
minimum drag endpoint. For the maximum volume endpoint there are strong shock waves, not only on the 
wing, as shown in Fig. 12, but also on the fuselage near the nose and tail. The fuselage shock waves are 
caused by the MOGA’s selection of minimum values for the x N and x T genes. These values help to 
maximize the vehicle volume, but cause shock waves to form at both the fuselage nose and tail because 
of excessive over-expansion. 

As the solution propagates along the Pareto front, trading volume to obtain a reduction in drag, the most 
noticeable initial change is associated with the fuselage nose and tail fairings. By choosing larger values 
for x N and x T , the drag is significantly decreased without sacrificing a large amount of volume. This can 
be seen by examining the intermediate solution in Fig. 13 and noting its position on the Pareto front in Fig. 
12. Further reductions in vehicle volume to achieve improvements in drag are achieved via the complex 
task of area ruling the fuselage in the vicinity of the wing-fuselage juncture. This situation can be observed 
by looking at the minimum drag solution in Fig. 13. More on this, including the comparison of selected 
fuselage cross-sections will be presented subsequently. 




A) Max volume solution B) intermediate solution C) Min drag solution 

Fig. 13 Surface Mach number contours for the wing/body optimization problem presented in Fig. 
12 at three positions along the Pareto front. 
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Pressure distributions along the fuselage for the three solutions described in Figs. 12-13 are presented in 
Figs. 14. Results are plotted along three fuselage circumferential stations— the keel line (cp = -90°), the 
side (cp = 0°), and the crown line (cp = 90°). Each of these solutions involves a high-mounted wing 
configuration, which is chosen by the MOGA optimizer regardless of the solution’s position along the 
Pareto front. More will be presented on this aspect subsequently. Because of this the side fuselage 
pressure distribution (cp = 0°) is plotted along a continuous circumferential line that lies below the wing for 
each of these three solutions. In addition, Mach number contours, as viewed from the fuselage’s side, are 
also displayed in Figs. 14. 



a) Maximum volume solution. 


b) Intermediate solution. 


1.5 



AXIAL DISTANCE ALONG FUSELAGE - X/L 

c) Minimum drag solution. 


Fig. 14 Fuselage pressure distributions for the wing/body optimization problem presented in Fig. 
12 at three positions along the Pareto front. 
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As shown in Fig. 14(a)— and previously discussed in conjunction with Fig. 13— the expansions, especially 
at the fuselage nose and tail, are quite large and lead to strong shock waves. As the fairings at the nose 
and tail are improved, see Figs. 14(b) and 14(c), the expansions and resulting shocks are greatly 
reduced. The shock wave induced on the fuselage by wing-fuselage interference can also be seen in Fig. 
14(a). Thinning of the wing, as well as the introduction of wing aft camber, reduces this effect in the 
intermediate solution [Fig. 14(b)]. Fuselage area ruling reduces the wing-fuselage interference shock 
strength even further as the minimum drag point is approached [Fig. 14(c)]. 

By comparing each part of Figs. 14 with the corresponding part in Figs. 13, it can be seen that the depth 
of the fuselage, as measured from keel to crown, is much larger than the width for any given longitudinal 
station. This is a direct result of the choices made in setting up the constraints for the A r array, as was 
already discussed in the fuselage parameterization section, and produces fuselage cross-sections that 
are elongated in the z-direction. This situation is presented in a more quantitative way in Fig. 15. 

A series of fuselage cross-sections are presented in Fig. 15 for the two-objective optimization problem 
displayed in Fig. 12. For each cross-section two curves are plotted. The solid curve corresponds to the 
solution obtained at the minimum-drag Pareto front endpoint, and the dashed curve corresponds to the 
solution obtained at the maximum-volume Pareto front endpoint. For this solution the wing root leading 
edge is placed at (x R , y R , z R ) =(1 .5, 0.0, 0.1 5) . Note that x R =1.5 and y R = 0.0 are values set at the 
beginning of the computation and are held fixed, and that z R = 0.15 is a value selected by the MOGA 
optimization for this particular chromosome from the range -0. 1 5 < z R < 0.1 5 (see Table 4). With this 
placement, the wing is high-mounted and intersects the fuselage between the second and third stations 
displayed in Fig. 15. Note: the fuselage cross-sections are drawn before the wing intersection is included. 




y/c 


y/c 


y/c 


y/c 


y/c 


Fig. 15 Selected fuselage cross-sections for the wing/fuselage optimization problem of Fig. 12 
showing differences between two Pareto front endpoints: C L = 0.45, M^= 0.84, ISELECT = 4, 

ITRAN = 0, N a =2, N c =34, N G =43, p = 0.1, ^=0.2, p 2 =0.2 and P = (0.04,0.32,0.32,0.32). 
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From this set of plots it is easy to see how the fuselage was area ruled in the vicinity of the wing-fuselage 
juncture to obtain the final amount of drag reduction. Volume is subtracted along the upper portion of the 
fuselage near the wing-fuselage intersection, creating “pear-shaped” fuselage cross-sections in this 
location. Every cross-section from the maximum volume solution contains more area than the 
corresponding minimum drag cross-section, because each of the A r values is individually maximized for 
the maximum volume solution. 

Convergence efficiency for the two-objective wing-fuselage optimization problem just discussed— drag 
minimization and volume maximization (both at fixed lift) — is presented in Figs. 16 and 17. Figure 16 
shows the previously discussed area error norm plotted versus the number of function evaluations for 
three different convergence histories involving three different population sizes— 34, 66, and 130 
chromosomes, respectively. Each curve is averaged over two ISEED values in an attempt to eliminate 
some of the statistical variation. Each of the convergence histories is nearly the same, indicating that 
convergence is not a function of population size. This is compatible with the results presented in Holst 
and Pulliam 15 where a number of model problems were used to evaluate GA convergence for two- 
objective optimization. 

For this set of computations the number of processors used to perform the function evaluations is equal to 
the number of chromosomes requiring function evaluations. There are always two passthrough 
chromosomes. The remaining chromosomes— those needing modification— are arbitrarily divided into 
three groups of equal size— or as close to equal as possible. One group is modified using crossover, the 
second using perturbation mutation, and the third using global mutation. Since two passthrough 
chromosomes are utilized for each population, the number of processors required to perform function 
evaluations is 32, 64, and 128, respectively. Using an increased number of chromosomes to define a 
given population allows an increased number of processors to be used for the optimization problem and 
thus, improves turn-around efficiency. 



NUMBER OF FUNCTION EVALUATIONS 


Fig. 16 Pareto front convergence history 
comparisons for three population sizes, each 
averaged over two ISEED values: C L = 0.45, 
M„- 0.84, ISELECT = 4, ITRAN = 0, N 0 = 2, 
N G =43, /? = 0.1, P! = 0.2, p 2 = 0.2. 



Fig. 17 Pareto front comparisons for three 
population sizes, each with two different 
ISEED values: C L = 0.45, M x = 0.84, 

ISELECT = 4, ITRAN = 0, N 0 =2, N G =43, 

13 = 0.1, p, =0.2, p 2 =0.2. 
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Figure 17 presents a comparison of the six Pareto fronts computed from the convergence efficiency study 
presented in Fig. 16. These results are statistically similar to previous results presented in this section and 
suggest that there is no significant effect of population size on the final solution. Care should be taken in 
evaluating this conclusion, as the number of solutions needed to construct a proper statistical average 
has not been utilized in this case because of computer time limitations. 

It is of interest to study how individual genes vary along the Pareto front. Selected genes from the wing- 
fuselage optimization problem described in Fig. 12 are plotted along the Pareto front in Figs. 18-19. They 
are plotted as a function of the Pareto front’s first objective, lift-to-drag ratio. Each gene is normalized so 
that it varies from 0 to 1 . For many genes the discrete points are plotted with a smoothed curve. In some 
cases, where the number of results on a single set of axes is large, only the smoothed curve is plotted. 



a) x -location of airfoil maximum thickness 
on lower surface at wing tip, GENE ., ^ =xlo 2 . 



c) Airfoil twist at wing root and tip, 

GENE 15 = Oi and GENE 16 = d 2 ■ 




b) Airfoil trailing edge angle at wing root and d) Vertical wing position (GENE 17 =z R ), 

tip, GENE 7 =ate., and GENE 14 = ate 2 . fuselage nose length (GENE 18 =x w ) and 

fuselage boattail length (GENE ig = x r ). 


Fig. 18 Normalized gene distributions along Pareto front from selected positions within each 
chromosome, wing/fuselage optimization problem: C L = 0.45, = 0.84, ISELECT = 4, ITRAN = 0, 

N 0 =2, A/ c = 34 , N g =43. 
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The first general observation is that a significant amount of statistical variation exists in many of the genes 
as they vary across the Pareto front. Some genes wildly oscillate without any discernible pattern. The x- 
location of the wing-tip lower-surface maximum thickness— xlo 2 , gene number 11— is such an example. 
Of the 43 available genes, xlo 2 contains the most variation. This gene is plotted in Fig. 17(a) and is seen 
to favor the maximum value in its allowable range. Nevertheless, it takes on many other values in a 
seemingly haphazard fashion. 

The reason for this behavior is straightforward. This gene has little impact on either objective being 
optimized and thus has little impact on the Pareto front. Any value is generally acceptable at any location 
along the Pareto front. To test this observation, several chromosomes were selected from the Pareto front 
for which xlo 2 was against one of its constraints. The value of xlo 2 within that chromosome was then 
moved to the opposite constraint without changing any other gene value, and a new solution was 
computed. The resulting design-space point moved little and was still a member of the Pareto front. 
Ideally, such a parameter shouldn’t be utilized as a gene in the optimization process. In reality it is difficult 
to know a priori which genes will behave in this way and which will not. 

Other parameters plotted in Fig. 18 include the trailing edge angle at the wing root and tip [Fig. 18(b)], the 
airfoil twist at the wing root and tip [Fig. 18(c)], and the vertical wing mounting position, the nose length, 
and the boattail length [Fig. 18(d)]. All of these gene parameters have reasonably controlled variations 
across the Pareto front. 

Figure 19 presents the variation of each of the A r parameters for the wing-fuselage optimization 
problem — all 24 of them — across the Pareto front. As already discussed, the A r parameters are 
organized into four circumferential stations, each with six longitudinally distributed values. Each 
circumferential station is plotted on a separate set of axes in Fig. 19. Figures 19(a), 19(b), 19(c) and 19(d) 
show the A r longitudinal variation along cp = -90°, -30°, 30° and 90°, respectively. 

As expected, all A r values at the maximum volume endpoint on the Pareto front lie at their maximum 
constraint limits. Moving away from this endpoint, many of the A r values decrease from their maximum 
constraint values and attain intermediate values. This is especially true for the A r values along the 
cp = 30° and cp =90° meridinal stations, where reductions in the various A r values have been determined 
by the MOGA optimizer to achieve optimal fuselage area ruling in the vicinity of the wing-fuselage 
juncture. 

It is of interest to study the effect of wing placement on the overall wing-fuselage optimization problem. In 
the cases presented above, the vertical wing position, z R , is constrained to lie between -0.15 and 0.15, 
while x R and y R are fixed at 1 .5 and 0.0, respectively. As a result of the optimization, the value of z R is 
chosen to lie at or near its upper constraint value, z R =0.15, over the entire Pareto front [see the 
variation of gene 17 plotted along the Pareto front in Fig. 18(d)]. In other words, a high wing intersection is 
favored by the optimization regardless of position on the Pareto front— at least in the context of the 
present design space construction. 

Figures 20 and 21 present results for a wing-fuselage optimization in which the wing intersection is 
constrained to be low mounted, that is, the value of z R is constrained to be -0.15. All other design space 
attributes and constraints are the same as in the previous wing-fuselage computations. This is 
implemented by setting the maximum constraint for z R equal to the minimum constraint, z R =-0.15, and 
by setting the corresponding z R masking array value to zero. Thus, for this problem, N G =42, instead of 
the previous value of 43. Note also that N c =66 for this series of computations, that is, 64 processors 
were used in the optimization. 

Figure 20 displays two computed Pareto fronts for the original wing-fuselage optimization problem and 
two Pareto fronts for the optimization that constrains the wing to be low-mounted. The two curves plotted 
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for each case correspond to two different ISEED values. As can be seen, the low mounted wing results 
are less efficient in the vicinity of the minimum-drag Pareto front endpoint. The low-mounted case 
achieves essentially the same level of drag, but only when vehicle volume is reduced. Note also that this 
effect is large in comparison to the statistical variation that exists between the two ISEED computations 
from each solution category. The effect of z R on the maximum volume Pareto front endpoint is negligible. 
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Fig. 19 Normalized gene distributions along the Pareto front from the fuselage portion of the 
parameterization: C L =0.45, M CO =0.84, ISELECT = 4, ITRAN = 0, N 0 = 2, A/ c =34, N G =43, /3 = 0.1, 
p A = 0.2 , p 2 = 0.2 , and P = (0.02, 0.33, 0.33, 0.32) . 


Also displayed in Fig. 20 are two circular symbols showing discrete chromosomes from each different 
category of Pareto front that most closely match a Lift-to-drag ratio of 13. Selected fuselage cross 
sections taken from these two points are compared in Fig. 21. A total of five longitudinal stations are 
presented with x/c = 0.49,1 .54,2.59, 3.64 and 4.69. The solid curve corresponds to the original case in 
which the wing vertical mounting position was allowed to float. (Although for this chromosome the value of 
z R that the MOGA selected was 0.15— a high mounted wing). The dashed curve corresponds to the case 
in which the wing was constrained to be low mounted. Despite being at the same L/D location on the 


34 


Pareto front, the two area rulings are quite different— quite different at each longitudinal station. Close 
examination shows how the MOGA modified the area ruling to account for the shifted wing mounting 
position, taking volume from the upper fuselage for the high-mounted case and from the lower fuselage 
for the low-mounted case. It is also obvious that the fuselage associated with the low-mounted wing 
contains a smaller amount of volume as is consistent with Fig. 20. 



Fig. 20 Pareto front comparisons for two different wing-fuselage optimization cases: a) baseline 
in which wing position is allowed to float (-0.15 <z R <0.15) and b) wing fixed in the low-wing 
position ( z R =-0.15). Each computation performed with two ISEED values: C L = 0.45, M ao = 0.84, 
ISELECT = 4, ITRAN = 0, N a =2, /V c =66, N G =43/42, 0-0.1, ^=0.2, p 2 = 0.2, and 

P = (0.02, 0.33, 0.33, 0.32). 
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Fig. 21 Selected fuselage cross-sections for the two wing/fuselage optimization problems of Fig. 
20 showing differences between the two solutions at LI D = 13: a) baseline in which wing position 
is allowed to float ( -0.15 < z R < 0.15) and b) wing fixed in the low-wing position (z R = -0.15): 
C L = 0.45 , = 0.84, ISELECT = 4, ITRAN-0, N 0 =2, /V c =66, N G =43/42, 0 = 0.1, p^O.2, 

p 2 = 0.2 , and P = (0.02, 0.33, 0.33, 0.32) . 
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Conclusions 


A multi-objective genetic algorithm (MOGA) optimization procedure, with several selection algorithm 
options and a new gene-space transformation procedure, has been presented and evaluated. It uses real- 
number encoding to represent all design space decision variables as genes and populations of fixed size 
to go from generation to generation. Four modification operators are utilized to advance from one 
generation to the next. They include passthrough, random average crossover, perturbation mutation, and 
mutation. The standard output for this approach is a Pareto front, which includes the best solutions from 
each objective, as well as a range of optimal “tradeoff” solutions in between. 

The MOGA optimization procedure was utilized to determine a number of multi-objective optimal solutions 
for two problems. The first involved a transonic wing lift-to-drag maximization (drag minimization) coupled 
with the simultaneous minimization of structural mass (all at fixed lift). The second involved a transonic 
wing-fuselage lift-to-drag maximization (drag minimization) coupled with the simultaneous maximization of 
vehicle volume (all at fixed lift). In every case the GA optimization process was convergent. However, 
some of the GA variations studied converged more rapidly than others. 

Of the two selection algorithms compared— greedy selection and bin selection— the latter produced the 
most efficient convergence across all of the problems studied. The gene-space transformation procedure 
produced mixed results in the area of GA convergence efficiency. In some cases it was moderately faster, 
for other cases moderately slower. 

The concept of aggressive mutation, which effectively biases the gene modification process toward gene 
values that reside on the design space boundary, was introduced. For the present problem, aggressive 
mutation produced a significant acceleration in convergence for the Pareto front. 

The new MOGA utilizes a masking array, which allows any gene (or set of genes) to be eliminated from 
the optimization process as a modifiable decision variable. This allows determination of the effect of a 
single gene (or gene subset) on the resultant Pareto front. 

A limited parallelization efficiency study was performed involving the use of 32 to 128 processors. 
Population sizes ranging from 34 to 130 chromosomes were used. The present MOGA is essentially 
embarrassingly parallel with regard to the number of chromosomes used for each generation, and the 
number of chromosomes did not affect convergence efficiency of the MOGA scheme. Thus, the 
turnaround time dramatically decreased with increasing number of processors. 
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